Tarea 5: MCMC

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Primera Parte: Preguntas Teóricas
  5. Segunda Parte: Elaboración de Código

Objetivo

a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

Pregunta 1: Model Evaluation and Information Criteria

Explique cómo cross-validation, criterios de información y regularización ayudan a mitigar o medir los problemas de underfitting y overfitting.

Respuesta

1. Cross-Validation:
Se dividen los datos en un conjunto de entrenamiento y validación para evaluar el desempeño del modelo con datos no vistos. Esto para: Detectar overfitting, ya que si un modelo que se ajusta demasiado al conjunto de entrenamiento tendrá un mal desempeño en los datos de validación.
Tambien identifica underfitting si el modelo no captura los patrones en ninguno de los conjuntos.

2. Criterios de Información :
Como el Akaike Information Criterion, evalúan el balance entre el ajuste del modelo y su complejidad: - Penalizan la complejidad excesiva (número de parámetros), ayudando a evitar overfitting. - Favorecen modelos con un buen ajuste pero simples, lo que reduce el riesgo de underfitting.

Por ejemplo, un modelo con menor AIC es preferido porque logra un equilibrio óptimo entre ajuste y parsimonia.

3. Regularización:
Incluye términos de penalización en la función objetivo del modelo, como en Ridge Regression. Esto controla los coeficientes del modelo:
- Reduce overfitting al penalizar coeficientes grandes, que podrían ajustarse demasiado a ruido en los datos.
- Evita underfitting permitiendo una flexibilidad adecuada para capturar patrones relevantes, especialmente si se optimizan hiperparámetros.

Pregunta 2: Directed Graphical Models

Diseñe una DAG para un problema causal inventado por usted de al menos 5 nodos (puede basarse en el ejemplo de Eugene Charniak) usando Dagitty y considere que la DAG tenga al menos: una chain, un fork y un collider. Muestre con dagitty todas las independencias condicionales de su DAG. Explique las independencias usando las reglas de d-separación.

Respuesta

# Cargar la librería
library(dagitty)
library(ggdag)
## Warning: package 'ggdag' was built under R version 4.4.2
## 
## Adjuntando el paquete: 'ggdag'
## The following object is masked from 'package:stats':
## 
##     filter
library(ggplot2)

# Crear la DAG
dag <- dagitty("
dag {
  Inteligencia -> Educacion -> Ingresos
  Inteligencia -> Oportunidades -> Ingresos
  Esfuerzo -> Ingresos
}
")

ggdag::ggdag(dag) + 
  ggtitle("DAG") +
  theme_dark()

  1. Educación ⊥ Esfuerzo | Ingresos:
    • Ingresos es un collider en el camino Educación -> Ingresos <- Esfuerzo.
    • Condicionar en Ingresos bloquea la dependencia entre Educación y Esfuerzo.
  2. Oportunidades ⊥ Educación | Inteligencia:
    • Inteligencia es un nodo de fork en el camino Educación <- Inteligencia -> Oportunidades.
    • Condicionar en Inteligencia bloquea el flujo de información entre Oportunidades y Educación.
  3. Oportunidades ⊥ Esfuerzo | Ingresos, Inteligencia:
    • Condicionar en Inteligencia bloquea la influencia desde Oportunidades a Ingresos.
    • Además, Ingresos como collider bloquea el camino hacia Esfuerzo.

Segunda Parte: Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
library(dplyr)
library(tidyr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)

# Manipulación de varios plots en una imagen.
library(gridExtra)

# Análisis bayesiano
library("rethinking")

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Pregunta 1: MCMC

El objetivo de esta pregunta es lograr samplear, mediante la técnica de MCMC, la distribución gamma.

En general la distribución gamma se denota por \(\Gamma(\alpha,\beta)\) donde \(\alpha\) y \(\beta\) son parámetros positivos, a \(\alpha\) se le suele llamar “shape” y a \(\beta\) rate La densidad no normalizada de una distribución gamma esta dada por:

\[ f(x\mid \alpha,\beta) = \begin{cases} x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\ 0 ~&\text{si } x \leq 0 \end{cases} \] donde \(\Gamma(\alpha)\) es una constante, usualmente se le llama función gamma.

    Escriba el algoritmo sin utilizar implementenaciones de la distribución gamma en r.

De ahora en adelante considere \(\alpha = 5\) y \(\beta = \frac{1}{5}\).

Respuesta

# Primeramente se define la función de densidad gamma

gamma_dens <- function(x, alpha, beta){
  if(x > 0){
    return(x^(alpha-1)*exp(-beta*x)) 
  }
  else{
    return(0)
  }
}

# Ahora se define el algoritmo de Metropolis Hasting con la función gamma anterior

metropolis_hasting <- function(theta_0, reps, alpha, beta){
  thetas <- vector(length = reps )
  current <- theta_0
  for ( i in 2:reps ) {
    thetas[i-1] <- current
    # distribución propuesta
    propuesta <- rnorm(1, mean = thetas[i-1], 1)
    # densidad gamma de la distribución propuesta
    gamma_propuesta <- gamma_dens(propuesta, alpha, beta)
    # densidad gamma de la distribución actual
    gamma_actual <- gamma_dens(thetas[i-1], alpha, beta)
    # Calculo del ratio
    ratio <- gamma_propuesta/gamma_actual
    # Probabilidad de aceptar
    prob <- min(ratio,1)
    decision <- rbinom(1,1,prob) 
    current <- ifelse( decision == 1 , propuesta , current )
  }
  return(thetas)
}

Luego, se consideran los valores de alpha y beta propuestos, y se tienen que \(\theta_{0} = 1\) variando la cantidad de \(reps\) quedará de la siguiente manera:

alpha <- 5
beta <- 1/5
theta_0 <- 1

# Variación de las reps
reps_list <- c(100, 1000, 10000, 1e5, 1e6)

par(mfrow = c(2, 2))

# Graficando
for (reps in reps_list) {
  muestras <- metropolis_hasting(theta_0, reps, alpha, beta)
  hist(muestras, probability = TRUE, main = paste("Repeticiones:", reps),
       xlab = "x", breaks = 50, col = "lightblue")
  
  #Para comparar con la distribución gamma real 
  curve(dgamma(x, shape = alpha, rate = beta), add = TRUE, col = "red", lwd = 2, lty = 1)

}

Es posible ver en los histogramas que inicialmente con una cantidad de repeticiones relativamente bajas no se logra apreciar la forma de una distribución gamma, sin embargo al llegar al orden de las 1e5 repiticiones ya se logra apreciar que es muy similiar y para 1e6 ya es prácticamente igual. Por lo que se puede concluir que para una mayor cantidad de repeticones el modelo de Metropolis-Hasting se va pareciendo cada vez más a una distibución gamma real.

Luego, para una cantidad de \(reps = 1e5\) y con variaciones del valor de \(\theta_{0}\). Se tendrá lo siguiente:

reps <- 1e5
# Variación de los thetas
thetas_list <- c(1, 2, 3, 4, 5, 6, 7, 8)

# Graficando
par(mfrow = c(2, 2))
for (theta in thetas_list) {
  muestras <- metropolis_hasting(theta, reps, alpha, beta)
  hist(muestras, probability = TRUE, main = paste("theta inicial:", theta),
       xlab = "x", breaks = 50, col = "lightblue")
  curve(dgamma(x, shape = alpha, rate = beta), add = TRUE, col = "red", lwd = 2, lty = 1)

}

Es posible ver que para todos los valores de \(\theta_{0}\) estudiados se tendrá que el modelo es bastante similiar al esperado, por lo que es posible concluir que está variable no tiene una gran influencia en el modelo de Metropolis-Hasting.

Finalmente se elegirán los valores \(\theta_{0} = 4\) y \(reps = 1e6\) por el análisis realizado anteriormente. Resultando:

theta_0 <- 4
reps <- 1e6

# Generando el modelo con Metropolis Hasting
muestras <- metropolis_hasting(theta_0, reps, alpha, beta)

# Graficar
hist(muestras, probability = TRUE, main = "Comparación con la distribución real",
     xlab = "x", breaks = 50, col = "lightblue")
curve(dgamma(x, shape = alpha, rate = beta), add = TRUE, col = "red", lwd = 2, lty = 1)
legend("topright", legend = c("Distribución  gamma real"), col = c("red"),
       lty = c(1, 2), lwd = 2)

Como es posible ver en el histograma, la aproximación obtenida con el modelo de Metropolis-Hasting coincide satisfactoriamente con la distribución gamma real. Por lo tanto este modelo daría como resultado estimaciones bastante acertadas con la realidad.

A work by CC6104